################################################################################
# REPLICATION SCRIPT FOR MAIN DOCUMENT
# "Far-right against Green: The re-emergence of geographically defined voting 
# patterns and the new environment cleavage in Western Europe"
# European Political Science Review
#
# Authors: Prof. Daphne Halikiopoulou, Dr. Christos Vrakopoulos, Dr. Christoph Arndt
# Origian Data: European Social Survey Round 8 (2016)
# We provide the cleaned data for the replication
################################################################################

# Clear workspace
rm(list = ls())

# Load required packages
required_packages <- c("nnet", "marginaleffects", "ggplot2", "gridExtra", 
                       "flextable", "dplyr", "officer")

for (pkg in required_packages) {
  if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
    install.packages(pkg)
    library(pkg, character.only = TRUE)
  }
}

################################################################################
# CUSTOM FUNCTIONS
################################################################################

# Two-way interaction plots
inter_plot <- function(
    model,
    variables,
    condition,
    type = "probs",
    groups = c("Far-Right", "Greens"),
    title = "",
    x_label = NULL,
    y_label = "Marginal effect",
    alpha = 0.1
) {
  if (is.null(x_label)) {
    x_label <- switch(condition,
                      "uerate_reg16" = "Regional Unemployment",
                      "gdpcap1000_reg16" = "Regional GDP",
                      condition)
  }
  
  pp <- plot_slopes(model, variables = variables, condition = condition, 
                    type = type, draw = FALSE)
  pp_filtered <- pp[pp$group %in% groups,]
  
  plot <- ggplot(pp_filtered, aes(x = .data[[condition]], y = estimate, group = group)) +
    geom_line() +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = alpha) +
    geom_hline(yintercept = 0, color = "black") +
    facet_wrap(~group, scales = "free_y") +
    labs(x = x_label, y = y_label, title = title) +
    theme_minimal() +
    theme(
      axis.title.x = element_text(size = 12),
      axis.title.y = element_text(size = 12),
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      legend.text = element_text(size = 12),
      plot.title = element_text(face = "bold", size = 12),
      strip.text = element_text(size = 12),
      legend.position = "none"
    )
  
  return(plot)
}

# Three-way interaction plots
inter_plot3 <- function(
    model,
    variables,
    condition,
    type = "probs",
    groups = c("Far-Right", "Greens"),
    title = "",
    x_label = NULL,
    y_label = "Marginal effect",
    alpha = 0.1
) {
  if (is.null(x_label)) {
    x_label <- switch(condition[1],
                      "uerate_reg16" = "Regional Unemployment",
                      "gdpcap1000_reg16" = "Regional GDP",
                      condition[1])
  }
  
  pp <- plot_slopes(model, variables = variables, condition = condition, 
                    type = type, draw = FALSE)
  pp_filtered <- pp[pp$group %in% groups,]
  
  plot <- ggplot(pp_filtered, aes(x = .data[[condition[1]]], y = estimate, 
                                  group = group, linetype = group)) +
    geom_line() +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = alpha) +
    geom_hline(yintercept = 0, color = "black") +
    scale_linetype_manual(values = c("solid", "dashed"), name = "Group") +
    facet_wrap(as.formula(paste("~", condition[2])), scales = "free_y") +
    labs(x = x_label, y = y_label, title = title) +
    theme_minimal() +
    theme(
      axis.title.x = element_text(size = 12),
      axis.title.y = element_text(size = 12),
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      legend.text = element_text(size = 12),
      plot.title = element_text(face = "bold", size = 12),
      strip.text = element_text(size = 12),
      legend.position = "right"
    )
  
  if (condition[1] == "uerate_reg16") {
    plot <- plot + scale_x_continuous(breaks = seq(2.5, 17, by = 2))
  }
  
  return(plot)
}

################################################################################
# LOAD DATA
################################################################################

# Load your cleaned ESS data here
# Example: ess <- readRDS("data/ess_cleaned.rds")

################################################################################
# DEFINE MODEL FORMULAS
################################################################################

# Direct effects models (Figure 1 and Table B1)
direct_formula1 <- party_vote ~ scale(energy_exp) + scale(uerate_reg16) + 
  domicil_4cat_fct + scale(hinctnta) + educ_level + Oesch8_fct + 
  gender + scale(age) + scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

direct_formula2 <- party_vote ~ scale(tax_ff) + scale(uerate_reg16) + 
  domicil_4cat_fct + scale(hinctnta) + educ_level + Oesch8_fct + 
  gender + scale(age) + scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

direct_formula3 <- party_vote ~ scale(coal_es) + scale(uerate_reg16) + 
  domicil_4cat_fct + scale(hinctnta) + educ_level + Oesch8_fct + 
  gender + scale(age) + scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

# Two-way interaction models (Figure 2 and Table B2)
h1_formula1 <- party_vote ~ scale(energy_exp) * scale(uerate_reg16) + 
  domicil_4cat_fct + Oesch8_fct + scale(hinctnta) + educ_level + 
  gender + scale(age) + scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

h1_formula2 <- party_vote ~ scale(tax_ff) * scale(uerate_reg16) + 
  domicil_4cat_fct + Oesch8_fct + scale(hinctnta) + educ_level + 
  gender + scale(age) + scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

h1_formula3 <- party_vote ~ scale(coal_es) * scale(uerate_reg16) + 
  domicil_4cat_fct + Oesch8_fct + scale(hinctnta) + educ_level + 
  gender + scale(age) + scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

# Three-way interaction models (Figure 3 and Table B3)
h2_formula1 <- party_vote ~ scale(energy_exp) * scale(uerate_reg16) * urbanrural + 
  Oesch8_fct + scale(hinctnta) + educ_level + gender + scale(age) + 
  scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

h2_formula2 <- party_vote ~ scale(tax_ff) * scale(uerate_reg16) * urbanrural + 
  Oesch8_fct + scale(hinctnta) + educ_level + gender + scale(age) + 
  scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

h2_formula3 <- party_vote ~ scale(coal_es) * scale(uerate_reg16) * urbanrural + 
  Oesch8_fct + scale(hinctnta) + educ_level + gender + scale(age) + 
  scale(lrscale) + scale(p_immi_index) +
  scale(trstplt) + scale(gdpcap1000_reg16) + scale(popden_reg16) + 
  scale(netmigr_reg16) + scale(tcoal_100jobs)

################################################################################
# PREPARE DATA
################################################################################

# Filter data
ess_analysis <- ess %>% 
  filter(countries == 1, Country != 20)

# Define color and shape values
color_values <- c("Far-Right" = "black", "Greens" = "darkgreen")
shape_values <- c("Far-Right" = 16, "Greens" = 17)

# Create output directory
dir.create("outputs", showWarnings = FALSE)

################################################################################
# FIGURE 1: DIRECT EFFECTS
################################################################################

# Model 1: Energy expenses
model_direct1 <- multinom(direct_formula1, 
                          data = ess_analysis, 
                          weights = anweight,
                          na.action = na.omit,
                          Hess = TRUE, 
                          model = TRUE,
                          trace = FALSE)

mb1_pred_ene <- marginaleffects::avg_predictions(model_direct1, by = "energy_exp")
mb1_pred_ene_fr_gr <- mb1_pred_ene %>% filter(group %in% c("Far-Right", "Greens"))

plot_mb1_pred_ene <- ggplot(mb1_pred_ene_fr_gr, 
                            aes(x = energy_exp, y = estimate, color = group, shape = group)) + 
  geom_point(position = position_dodge(0.25)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1, 
                position = position_dodge(0.25)) +
  scale_color_manual(values = color_values) +
  scale_shape_manual(values = shape_values) +
  scale_x_continuous(
    breaks = 1:5,
    labels = c("Not at all worried", "Not very worried", "Somewhat worried", 
               "Very worried", "Extremely worried")
  ) +
  labs(title = "Energy becomes too expensive",
       x = "", y = "Predicted Probability",
       color = "Party Family", shape = "Party Family") +
  theme_minimal() +
  theme(legend.position = "bottom",
        axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        plot.title = element_text(face = "bold", size = 12)) +  
  coord_flip()

# Model 2: Fossil fuel taxes
model_direct2 <- multinom(direct_formula2, 
                          data = ess_analysis, 
                          weights = anweight,
                          na.action = na.omit,
                          Hess = TRUE, 
                          model = TRUE,
                          trace = FALSE)

mb2_pred_tax <- marginaleffects::avg_predictions(model_direct2, by = "tax_ff")
mb2_pred_fr_gr_tax <- mb2_pred_tax %>% filter(group %in% c("Far-Right", "Greens"))

plot_mb2_pred_tax <- ggplot(mb2_pred_fr_gr_tax, 
                            aes(x = tax_ff, y = estimate, color = group, shape = group)) + 
  geom_point(position = position_dodge(0.25)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1, 
                position = position_dodge(0.25)) +
  scale_color_manual(values = color_values) +
  scale_shape_manual(values = shape_values) +
  scale_x_continuous(
    breaks = 1:5,
    labels = c("Strongly in favour", "Somewhat in favour", "Neither in favour nor against", 
               "Somewhat against", "Strongly against")
  ) +
  labs(title = "Support for higher taxes on fossil fuels",
       x = "", y = "Predicted Probability",
       color = "Party Family", shape = "Party Family") +
  theme_minimal() +
  theme(legend.position = "bottom",
        axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        plot.title = element_text(face = "bold", size = 12)) +  
  coord_flip()

# Model 3: Coal support
model_direct3 <- multinom(direct_formula3, 
                          data = ess_analysis, 
                          weights = anweight,
                          na.action = na.omit,
                          Hess = TRUE, 
                          model = TRUE,
                          trace = FALSE)

mb3_pred_coal <- marginaleffects::avg_predictions(model_direct3, by = "coal_es")
mb3_pred_fr_gr_coal <- mb3_pred_coal %>% filter(group %in% c("Far-Right", "Greens"))

plot_mb3_pred_coal <- ggplot(mb3_pred_fr_gr_coal, 
                             aes(x = coal_es, y = estimate, color = group, shape = group)) + 
  geom_point(position = position_dodge(0.25)) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1, 
                position = position_dodge(0.25)) +
  scale_color_manual(values = color_values) +
  scale_shape_manual(values = shape_values) +
  scale_x_continuous(
    breaks = 1:5,
    labels = c("Strongly in favour", "Somewhat in favour", "Neither in favour nor against", 
               "Somewhat against", "Strongly against")
  ) +
  labs(title = "Support for coal as energy source",
       x = "", y = "Predicted Probability",
       color = "Party Family", shape = "Party Family") +
  theme_minimal() +
  theme(legend.position = "bottom",
        axis.title.x = element_text(size = 12),
        axis.title.y = element_text(size = 12),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        plot.title = element_text(face = "bold", size = 12)) +  
  coord_flip()

# Combine and save Figure 1
figure1 <- grid.arrange(plot_mb1_pred_ene, plot_mb2_pred_tax, plot_mb3_pred_coal, nrow = 3)
ggsave("outputs/Figure1_direct_effects.png", figure1, 
       width = 12, height = 12, units = "in", dpi = 300)

################################################################################
# FIGURE 2: TWO-WAY INTERACTIONS
################################################################################

# Model H1.1: Energy × Unemployment
model_h1_1 <- multinom(h1_formula1, 
                       data = ess_analysis, 
                       weights = anweight,
                       na.action = na.omit,
                       Hess = TRUE, 
                       model = TRUE,
                       trace = FALSE)

plot_h1_1 <- inter_plot(
  model_h1_1,
  variables = "energy_exp",
  condition = "uerate_reg16",
  x_label = "Regional Unemployment",
  title = "Energy becomes too expensive"
)

# Model H1.2: Taxes × Unemployment
model_h1_2 <- multinom(h1_formula2, 
                       data = ess_analysis, 
                       weights = anweight,
                       na.action = na.omit,
                       Hess = TRUE, 
                       model = TRUE,
                       trace = FALSE)

plot_h1_2 <- inter_plot(
  model_h1_2,
  variables = "tax_ff",
  condition = "uerate_reg16",
  x_label = "Regional Unemployment",
  title = "Support for taxes on fossil fuels"
)

# Model H1.3: Coal × Unemployment
model_h1_3 <- multinom(h1_formula3, 
                       data = ess_analysis, 
                       weights = anweight,
                       na.action = na.omit,
                       Hess = TRUE, 
                       model = TRUE,
                       trace = FALSE)

plot_h1_3 <- inter_plot(
  model_h1_3,
  variables = "coal_es",
  condition = "uerate_reg16",
  x_label = "Regional Unemployment",
  title = "Support for coal as energy source"
)

# Combine and save Figure 2
figure2 <- grid.arrange(plot_h1_1, plot_h1_2, plot_h1_3, nrow = 3)
ggsave("outputs/Figure2_twoway_interactions.png", figure2, 
       width = 12, height = 12, units = "in", dpi = 300)

################################################################################
# FIGURE 3: THREE-WAY INTERACTIONS
################################################################################

# Model H2.1: Energy × Unemployment × Place
model_h2_1 <- multinom(h2_formula1, 
                       data = ess_analysis, 
                       weights = anweight,
                       na.action = na.omit,
                       Hess = TRUE, 
                       model = TRUE,
                       trace = FALSE)

plot_h2_1 <- inter_plot3(
  model = model_h2_1,
  variables = "energy_exp",
  condition = c("uerate_reg16", "urbanrural"),
  title = "Energy becomes too expensive"
)

# Model H2.2: Taxes × Unemployment × Place
model_h2_2 <- multinom(h2_formula2, 
                       data = ess_analysis, 
                       weights = anweight,
                       na.action = na.omit,
                       Hess = TRUE, 
                       model = TRUE,
                       trace = FALSE)

plot_h2_2 <- inter_plot3(
  model = model_h2_2,
  variables = "tax_ff",
  condition = c("uerate_reg16", "urbanrural"),
  title = "Support for taxes on fossil fuels"
)

# Model H2.3: Coal × Unemployment × Place
model_h2_3 <- multinom(h2_formula3, 
                       data = ess_analysis, 
                       weights = anweight,
                       na.action = na.omit,
                       Hess = TRUE, 
                       model = TRUE,
                       trace = FALSE)

plot_h2_3 <- inter_plot3(
  model = model_h2_3,
  variables = "coal_es",
  condition = c("uerate_reg16", "urbanrural"),
  title = "Support for coal as energy source"
)

# Combine and save Figure 3
figure3 <- grid.arrange(plot_h2_1, plot_h2_2, plot_h2_3, nrow = 3)
ggsave("outputs/Figure3_threeway_interactions.png", figure3, 
       width = 12, height = 12, units = "in", dpi = 300)

# Replication complete - all outputs saved to 'outputs/' directory